Dynamical Critical Behaviors of the Ising Spin 
Chain: Swendsen-Wang and Wolff Algorithms 

P. L. Krapivsky 

Center for Polymer Studies and Department of Physics, Boston University, Boston, 
MA 02215, USA 

E-mail: paulk@bu.edu 
Abstract. 

We study the zero-temperature Ising chain evolving according to the Swendsen- 
Wang dynamics. We determine analytically the domain length distribution and various 
"historical" characteristics, e.g., the density of unreacted domains is shown to scale 
with the average domain length as (l)~ s with 5 = 3/2 (for the g-state Potts model, 
5 = 1 + q^ 1 ). We also compute the domain length distribution for the Ising chain 
endowed with the zero-temperature Wolff dynamics. 



PACS numbers: 02.50Ey, 05.40.+j, 82.20.Mj 
1. Introduction 

Interesting and relatively poorly understood dynamical critical behaviors occur when 
statistical-mechanical systems are quenched from a disordered phase to their critical 
points. For the Ising spin systems, two very popular dynamics were introduced long ago 
by Glauber [T| and Kawasaki [2|. Glauber and Kawasaki algorithms are the simplest 
dynamical rules based on local moves — single spin flips for the non-conservative Glauber 
dynamics and spin exchanges for the conservative Kawasaki dynamics. Glauber's and 
other single spin-flip dynamics, particularly the Metropolis algorithm, have become a 
powerful tool for understanding the equilibrium behavior of the statistical-mechanical 
systems well away from the critical temperature [Hill]. The simulation becomes very 
slow at criticality, however, and to overcome this difficulty non-local moves, such as 
cluster flips, have been suggested. The Swendsen-Wang [5] and Wolff [0] dynamics are 
two well-known cluster algorithms that are widely used in elucidating the equilibrium 
critical behavior in statistical physics and lattice field theory. The dynamical critical 
behavior of these algorithms is an active area of research, see [7] and references therein. 
This work heavily relies on simulations, e.g., the value of the dynamical critical exponent 
for the two-dimensional Ising model endowed with the Swendsen-Wang dynamics is 
known only numerically; analytic studies of the Swendsen-Wang dynamics have been 
limited so far to the Ising model on the complete graph, that is to the Curie- Weiss or 
the mean- field model El UD] • 
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The purpose of this paper is to investigate the dynamical aspects of the cluster 
algorithms — particularly the Swendsen- Wang algorithm — in the simplest possible 
setting, that is in one dimension. The critical temperature is usually T c = for one- 
dimensional systems. The zero-temperature dynamics can be quite peculiar, for example 
the Ising spin chain subject to the zero-temperature Kawasaki dynamics freezes [TT|I12|. 
Non-conservative dynamics, however, usually bring the system to the ground state, e.g., 
the zero-temperature Ising-Glauber chain reaches the ground state in a time r ~ C 2 
(here £ is the system size), that is the dynamical exponent is z = 2 for the Glauber 
algorithm QjJ [I^j. We will see that z = for the Swendsen- Wang algorithm in one 
dimension; more precisely, r ~ ln£. The one-dimensional version of the Swendsen- Wang 
algorithm is also an appealingly simple model that is reminiscent of other important 
models of phase-ordering dynamics like the time-dependent Ginzburg-Landau equation 
with no thermal noise (i.e., T = 0) f3JEI|. Furthermore, the one-dimensional Swendsen- 
Wang algorithm provides a useful laboratory to probe not merely the dynamical critical 
exponent but various much more subtle dynamical characteristics. 

This paper is organized as follows. In section 2 we show that the Ising chain 
endowed with zero-temperature Swendsen- Wang dynamics exhibits scaling with the 
average length growing exponentially with time. Section 2 also contains the derivation 
of the domain length distribution and a number of subtle statistical properties of the 
domains like the density of domains which never flipped. In section 3 we investigate 
the Ising chain endowed with zero-temperature Wolff dynamics. This model has been 
previously studied by Derrida and Hakim [EJ; here we further analyze the model, 
particularly we determine the domain length distribution. A summary is given in the 
last section 4. 

2. Swendsen- Wang Dynamics 

In one dimension, the Ising spin chain can be thought as an array contiguous alternating 
domains of up and down spins. At zero temperature, the Swendsen- Wang dynamics 
amounts for randomly choosing a domain and flipping it. The flipping of a domain 
I implies that it merges with its left and right neighboring domains i_ and I + thus 
forming a single domain 7_ U I U I + . 

2.1. Domain Length Distribution 

Let Ni(t) is the number of domains of length I and N(t) = J2i>i N[(t) is the total number 
of domains. The average number of domains that flip in an infinitesimal time interval 
At is equal to N(t) At. In every flipping event three domains merge into a single one, 
so N(t) changes according to 

N(t + At) = N(t) - 2At N(t). (1) 
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Similarly, Ni(t) evolves according to 

N l (t + At) = N l (t)-3AtN l (t) + At ^ ^ Nj {t)^. (2) 

The term 3 At Ni(t) accounts for the loss that occurs when the domain or either of its 
neighbors is flipped, while the last term on the right-hand side of Eq. (j2J) accounts for 
the gain due to the flipping of a domain of length j followed by immediate merging 
with two adjacent domains of lengths i and k. Equation $1} is obviously exact. The 
linear loss term in (J2J) is also exact, while the non-linear gain term is exact only if 
the sizes of adjacent domains are uncorrelated. However, whenever a domain merges 
with two adjacent domains, the resulting domain does not acquire any correlation with 
the neighbors, i.e., correlations are not dynamically generated. (See Ref. jT3] for a 
detailed justification in the context of a somewhat similar model, viz. the noiseless 
time-dependent Ginzburg-Landau equation in which domains also merge with their 
neighbors.) Therefore, if initially the domain lengths were uncorrelated, they remain 
uncorrelated at all later times and Eq. (J2J) is exact. 

The total length of the system is C = J2i>ilNi(t). In the thermodynamic limit 
C — > oo, it is convenient to use the domain length densities rty(i) = Ni(t)/£ and the 
(total) domain density n{t) = J2i>i n i(f) — N(f)/C From ((Tj) we find that the domain 
density evolves according to 

dn „ . . 

Solving (J5J) gives n(t) = n(0) e~ 2t . Therefore the average domain size (/) = 1/n increases 
exponentially. In contrast, the average size exhibits an algebraic growth (I) ~ t 1 ^ [17 j 
in most models describing domain coarsening following a quench to zero temperature. 
Likewise, the equation for Ni(t) leads to 
dm ^ niUjUh 

The form of Eqs. (J3J) suggests to consider normalized densities pi(t) = ni(t)/n(t). These 
quantities satisfy 

~rr = ~Pi + E PiPiPk- ( 5 ) 
ai i+j+k=i 

Equations (HJ)-© are mathematically similar to equations describing the 3-particle 
coalescence process [THJ and can be solved accordingly. Introducing the generating 
function 

R{t,x) =J2x l pi(t) (6) 
i>i 

we convert an infinite system of equations (0) into a single equation 

£ = + (T) 
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Solving (J7J) gives 



fl(t,*)= , e ~ tR ° (x) (8) 



with Rq(x) = -R(0,x). For instance, consider the evolution starting from the highest 
energy antiferromagnetic state. In this case, pi(0) = 6n, i.e. Rq(x) = x. Inserting this 
into (JHJl and expanding in powers of x we obtain p2i{t) = and 

p 2 mW = e-(^f^y( 2 /). (9) 
The original densities read 

W0 = c- at ( iz f^)'( 2 /)- (io) 

In the scaling limit /, t — > oo with the scaling variable 

L = /e~ 2t = finite, (11) 
the densities become 

n / (t) = e- 4 ^(L), ^(L) = exp(-L/2). (12) 

For other initial conditions it is quite difficult to extract explicit results for pi from the 
general solution (jSj) for the generating function. The most natural situation arises when 
the system at T = oo is suddenly quenched to T = 0. The appropriate initial condition 
is completely uncorrelated with each spin taking the values ±1 independently and with 
equal probabilities. Then pi(0) = 2~ l , or Rq(x) — x/(2 — x) thus leading to a relatively 
simple exact expression for the generating function. The resulting pi admits a neat 
expression 



ft(t)=T . 2Fl ( / ,L 1; i^|), _i±yjZ£Z (13) 

An apparent simplicity of this solution is illusory as (JT3j) involves the hypergeometric 
function. Fortunately, the most interesting scaling behavior ()12j) is universal, i.e., 
independent on initial conditions (modulo the assumption that the decay p/(0) vs. I is 
sufficiently steep). Therefore in the following we focus on the antiferromagnetic initial 
condition. 

Note that the domain length distribution n/(t) is beyond the reach of analytical 
approaches [Hj for the majority of models of domain coarsening; e.g., for the Ising 
chain endowed with zero-temperature Glauber dynamics the distribution n\(t) is still 
unknown although a few exact and approximate results are available ^01 120] • The 
Swendsen- Wang dynamics is obviously more tractable than the Glauber dynamics - 
we determined ni(t) for the antiferromagnetic initial condition and an exact scaling 
expression in the general case. 
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2.2. Domain Number Distribution 

For the Swendsen- Wang dynamics it is also possible to probe analytically various 
historical characteristics. The simplest such quantity is the density m p (t) of domains 
composed of p 'parent' domains that never flipped. (Each such domain can of course 
include domains that flipped during the time interval (0, i).) The domain number 
distribution m p (t) is formally defined as follows [2D]- Initially m p (0) = 8 P \. In every 
merging event, the central domain flips while the two adjacent domains do not flip; 
therefore if they have a and j3 parents, respectively, the resulting domain has a + (3 
parents. The domain number distribution m p (t) evolves as follows: 
dm p x - m a mf3 

~d = + £ — • (14) 

A solution to these equations has an exponential form 

m p = Aa p ~ 1 . (15) 

This ansatz reduces an infinite system (JUjl to a couple of differential equations 

dA n . da i . ,.„ s 
M = - iA ' Hi = ' < 16 > 

where n(t) = e~ 2t , see Eq. (jSJ). Solving (fTBj) subject A(0) = l,a(0) = (implied by the 

initial condition m p (0) = S p i) we finally obtain 



m p (t)=e- 3t (l-e- t Y '. (17) 



p-i 



Thus P = pe * is the scaling variable and the scaling form is simply exponential: 

m p {t) = e~ 3t g(P), g(P) = e - p . (18) 

The survival of a domain with and without merging are characterized by the 
domain persistence exponents ip and 5 which were investigated for the Ising-Glauber 
|2D] and Ising-Kawasaki j2I] spin chains. (Domain persistence and closely related cluster 
persistence exponents were also studied for several other models [22 123 EH 123 ) The 
exponent 5 describes the decay of primordial domains 

m 1 ~(/)-' 5 . (19) 

The exponent if counts the average number of parents per domain (p) = YliV m vl S m p : 

<P> - (20) 
In the present case, m\ = e~ 3t and (p) = e*, while the average domain length is (I) = e 2t . 
Therefore 

t = I * = \ (2D 

Even for the simplest models, these exponents are still known only numerically, e.g., 
5 m 2.54 and ip w 0.252 for the Ising-Glauber spin chain |2D| . 5 w 2.12 and ^ w 0.39 for 
the Ising-Kawasaki spin chain [I^j. For a few models, however, the domain persistence 
exponents have been computed analytically [2011221123, e.g., for the random field Ising- 
Glauber spin chain 5 = oo and ip = (3 - V§)/4 = 0.190983 ... [231- 
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2. 3. Domain Length-Number Distribution 

The (normalized) domain length-number distribution pi p captures both the spatial 
and historical characteristics of the coarsening domain mosaic and contains previous 
distributions: p\ = J2 P Pip an d m p = nJ^iPip- The domain length-number distribution 
satisfies 

= -Pip + J2 PiaPjPkf3- (22) 

i+j+k=l a+j3=p 

To determine p\ v we use the two-variable generating function 

n(t,x,y) = j2J2 xl y p pip( t )- ( 23 ) 

1>1 P>1 

Multiplying Eq. f!22l) by x l y p , summing over all l,p > 1, and using already known result 
(JSJ) for the one- variable generating function J2j X ' J Pjif) we find that 7?. satisfies 

dt Ji _ (i _ e -2fW2 



Solving this equation subject to 7^.(0, x, y) = xy we obtain 

n(t lXt y) = e ~ txV (25) 

i -y + - (l -e- 2 ')x 2 

Expansion in y is simple and for every p we get 

p-i 



X; =e"'x e" 2 *)x 2 . (26) 

The expansion in x is also straightforward and it leads to a series representation for pi p . 
Of course, pi p = for even Vs. For odd Vs we find 

P2i-i, P =0 for / < p, 
P2 P -i,p — e , 



1 _ (I -e- 2t \ p 
p2 P +i,p = - (p - 1) e * I J , 

1 _ fl-e- 2t \ p+1 
P2 P +3,p= - (p- l)(p + 2)e * I j , 

etc. To extract the scaling behavior, it is more convenient to use (J2"H|) rather than ()2fij) . 
Writing 

x = 1 - e" 2 ^, y = 1 - e-\ (27) 
and taking the limit t — > oo we simplify (|25|) to 

n ^v) = —-^= W - (28) 
77 + VI + 2£ 

In the scaling limit l,p,t — > 00 with the scaling variables 

L = I e~ 2t = finite, P = pe~ t = finite, (29) 
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the domain length-number distribution admits the scaling form 

p lp (t) = e- 3t H(L,P). (30) 

Inserting (|2"7|). (|3*0)l into Eq. (|2"3*)) and replacing summation by integration we convert 
the two- variable generating function into the double Laplace transform 

K(£,V)= dLe^ L dP e~ vP H(L, P) . (31) 
Jo Jo 

Using Eqs. ()28j1 . (|31jl and performing the inverse Laplace transform we obtain 

H(L - P) ^ exp H-5l)- < 32 > 

Comparing (J3*2"j) with individual distributions (fT2|) and (|T%|) we see that the domain 
length- number distribution does not factorize even in the scaling limit. 

From the length-number distribution one can extract a lot of things, e.g. the 
fraction of persistent spins, i.e., spins which have not flipped during the time interval 
(0, t); this quantity usually decays as (l)~ e , where 9 is the persistence exponent (2E|- For 
the antiferromagnetic initial condition, for instance, the number of persistent spins in a 
domain is exactly equal to the number of parents. The average number of parents is 

<p>i = ^^ = pr x £p*P- (33) 

2^p>i Pip p >i 

In the long time limit, we can use (J32|) and replace the summation by integration. This 
leads to the asymptotic 

— when I — * oo. (34) 

The fraction of persistent spins n^ 1 J2i>i(p) i n i is now computed to give e~* = (/) _1 / 2 . 
Thus, the persistence exponent is 6 = 1/2. 

2-4- The q-state Potts model 

Some of the above calculations can be generalized to the case of the Potts model. For 
the g-state Potts model, each domain is in one of the q possible states, and each time a 
domain is updated it adopts the state of one of the two adjacent domains. The updating 
of a domain results in merging of all three domains with probability l/(q — l) while with 
probability (q — 2)/(g — 1) only two domains merge. Therefore the average number of 
domains lost in every merging event is 2 x + 2=j| = -2L leading to 

dn q 

— = n. (35) 

dt q - 1 ( J 

Therefore the average domain size (/) = 1/n increases exponentially for arbitrary q. 

The number distribution m p (t) is also readily computable for the g-state Potts 
model. One has 

dm p g + 1 1 m a m p 

—r~ — 7^»H 7 > -■ (36) 

dt g-1 p g-1 n v; 

^ ^ a+p=p 
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This equation admits the exponential ansatz ()15j) that reduces into a couple of 
differential equations 



dA 
~dt 



q + 1 



q-1 



A, — 



da 



(g — f )n ' 



where n = exp 



— M 

9-1 



Solving (J3TJ) subject to the initial conditions A(0) 



a(0) = (implied by m p (0) = <5 p i) we obtain 

m P (t) = Q 9+1 (i-Q) p - 1 

Re-expressing the quantities mi = Q q+1 and (p) 



g(t) = e -*/(«/-i). 



(37) 
1 and 

(38) 



through the average domain 
length (I) = n~ x = Q~ q we find that the domain persistence exponents defined via 
equations (fTT)|) - (f2T)|) are given by 

. q + 1 . q — 1 



(39) 



q q 

The exponents thus obey the sum rule 6~(q) + ip(q) = 2. 

Now consider the length distribution. Particularly, the normalized domain length 
distribution evolves according to 

dpi 1 V- L ?- 2 

■" u i i+j+ k=i y 1 i+j=i 



q 

From this equation, we find an implicit relation for the generating function 



(40) 



[l-Ry- 1 (R + q-l) 
Ri 



,9* 



R^- 1 (R + q 



Rq 



This relation is a polynomial of R of degree q. Hence it is impossible to determine an 
explicit relation for the generating function, R(t,x) = Rq(x)], apart from the Ising 
case (q = 2) and two next cases q — 3,4. The explicit expressions for q = 3, 4 are very 
involved so the exact computation of ni(i) looks daunting. 

Rather than seeking an exact solution, let's consider the asymptotic behavior. The 
scaling ansatz 



Pi(t) 
recasts (jjUj) into 



nl, 



J r + qL^ + J r *J r *J r 
dL 



{q-2)F*F=0 



(41) 



(42) 



where symbol * denotes the convolution operation, so that T*T is the usual convolution 
/ L dL x T{L X ) T{L - Li), and T * T * T = / L dL x dL 2 Ffa) T(L 2 )T(L -L x - L 2 ). 
The Laplace transform $(s) = / °° dLe~ sL T{L) satisfies 

d$ 



qs— = $ d + (q - 2)$ z 
ds 

whose (implicit) solution reads 

"$ + q - 1 



9-1 



(43) 



(44) 
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The sum rules J2 Pi — 1 an d Pi = n 1 l ea d to 

/ dLJ r (L) = l, / dLLT{V) — 1. (45) 
Jo Jo 

These two constraints determine first two constants in the small s expansion of the 

Laplace transform: $(s) = 1 — s + This behavior was taken into account in fixing 

a constant in the general solution to Eq. (|4*H|) . 

To complete the task, we must find from ()44|) and then perform the inverse 

Laplace transform. The first step therefore requires finding a root of the polynomial of 

$ of degree q. Thus it is apparently impossible to find an explicit scaling solution when 

q > 5. However, we can deduce the most interesting asymptotics for an arbitrary q. For 

instance, Eq. (|44jl yields $ — > (1 — q^ 1 ) 1 ^ s~ 1+1 / Q as s — > oo, from which 

■F(L) -> V ? gy , L-~> as LjO. (46) 

The large L behavior of T{L) is reflected by the type of the (closest to the origin) 
singularity of its Laplace transform. From Eq. (J44j) we find that the singularity is the 
branch point located at s* = — g -1 /^ -1 ), namely 

_1 q-2 r 1_ 

$^22 g2(,-i) s + q 9-1 
leading to 

■F(L)^g^T exp\-Lq~^] (47) 

as L —>■ oo. Thus the large L behavior is qualitatively similar for all q, while the small 
L behavior is affected by the number of states. 

The Swendsen- Wang dynamics of the g-state Potts model is particularly simple in 
the q — ► oo limit when only two adjacent domains can merge. The domain number 
distribution is trivial in this case, m p (t) = e _< 5 P \. The normalized domain length 
distribution satisfies 

^ = -Pi + E PiPr (48) 

i+j=l 

This system of equations resembles (JUj) and the solution is accordingly found by seeking 
pi(t) in an exponential form like (|15j). For initially uncorrelated Potts spins, all domains 
have initial length one when q = oo. Therefore pi(0) = 5n and the solution reads 

p,(*) = e" t (l-e" t ) , ~ 1 - ( 49 ) 
The scaling form of the domain length distribution is exponential, T{V) = e~ L ; of 
course, this result can also be derived from (J44j) which in the q — > oo limit simplifies to 
$ = (1 + s) _1 . An apparent contradiction of the pure exponential scaled domain length 
distribution and the general large L asymptotic (|4Tj) is the indication that the limits 
q — » oo and L — >• oo do not commute. 
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3. Wolff Dynamics 

At zero temperature, the Wolff dynamics amounts for randomly choosing a spin and 
flipping the whole domain containing that spin. The flipping of a domain / again 
implies that it merges with its left and right neighboring domains J_ and I + to form a 
domain I U / U I + . In contrast with the Swendsen- Wang dynamics, the flipping of a 
domain now occurs with rate proportional to its length. Therefore the domain density 
decreases with constant rate 

£ = < 50 » 

i.e., n{t) = 1 — It and the whole system reduces to a single domain at t c = 1/2. 
The governing equations for the domain length densities read 

— = -lni-2—+ ^ 1— . (51) 

We again use the normalized densities pi(t). They satisfy 

^ = -I pi + JPiPjPk- (52) 

at i+j+k=i 

The generating function (JOJ) satisfies 

— = -a:— + xir— . (53) 
or ax ax 

Changing variables from (t, x) to (T,X) = {t,t — Inx) removes the first term on the 

right-hand side of Eq. (J53*Jl : 

»5 = -J?** (54) 
Rewriting for X = T(R,X) gives 

= ^ (55) 
which is solved to yield X = F(R) + R 2 T, or 

t-lnx = F(R) + R\ (56) 

with F(R) depending on initial conditions. For the antiferromagnetic initial condition 
Ro(x) = x and thence F(R) = — \nR, so the exact solution (}56|) becomes 

x = Re t-R 2 t. (57) 

Clearly, Eq. (|57j) gives an expansion of x in terms of i?. We are seeking the opposite, 
R = J2 pix 1 . To determine pi we write 

which is combined with Eq. ()57|) to give 



M*) = dR 



lR?t plR^t 

-2t 



R l R l ~ 



(58) 
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An elementary computation shows that p2i = and 




(59) 



Since ni = npi and n — 1 — 2t, the densities vanish at t c — 1/2. Note also that in the 
proximity of the critical point (i.e., when n — > and / — > oo), Eq. ()59|) simplifies to 



This expression was previously derived in Ref. fB] by a direct analysis of the generating 
function near the critical point. Equation (p]j) shows that (/) = n~ x does not characterize 
the domain length distribution: Almost all domains are short with lengths of order 1 
but because the domain length distribution is algebraic, pi oc /~ 3 / 2 with a cutoff length 
of the order of n~ 2 , the average domain length diverges as n~ l . 

The single domain covers the entire spin chain at t c = 1/2, i.e., the system undergoes 
a gelation transition. This transition differs from the ordinary gelation transition that 
occurs in mean-field models of polymerization j23 |2Hj and evolving random graphs [29] 
where the giant component that is born at the critical time undergoes a long period 
of growth before it engulfs the entire the system. The difference from the ordinary 
gelation transition is not surprising: In one dimension, the giant component must cover 
the entire system, so the transition is discontinuous, while in the mean-field models 
the transition is continuous. Despite of this important physical distinction, the domain 
length distribution (joT^ is very similar to the cluster size distribution in polymerization 
[28j and the component size distribution in evolving random graphs [29J. 

4. Summary 

We demonstrated that the Ising chain endowed with zero-temperature Swendsen- Wang 
dynamics exhibits scaling with the average length growing exponentially with time. We 
computed the domain length distribution in special cases, e.g. for the antiferromagnetic 
initial condition and the random initial condition. The scaled domain length distribution 
was shown to be a product of a power-law and an exponential over the entire range of the 
(scaled) length. We also computed the domain number distribution, the domain length- 
number distribution, and the domain persistence exponents. The domain length-number 
distribution does not factorize into the product of the single variable distributions. Some 
of the calculations have been generalized to the Potts model; in particular, the domain 
number distribution and the domain persistence exponents have been obtained. 

We also studied the Ising chain endowed with zero-temperature Wolff dynamics 
and showed that the system undergoes a kind of gelation - below the critical time the 
total number of domains is an extensive variable (proportional to the system size) while 
at the critical time the entire system gets covered by the single domain. This gelation 
transition is discontinuous since in one dimension, the giant component must cover the 
entire system; in contrast, in the mean- field models the gelation transition is continuous. 
Still we found that the domain length distribution in the zero-temperature Ising- Wolff 
chain is mathematically similar to the cluster size distribution in polymerization. 




(60) 
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